Vortex nucleation in rotating Bose-Einstein condensates 
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We study the formation and stability of a single vortex state in a weakly-interacting Bose-Einstein 
condensate that is confined in a rotating harmonic potential. Our results are consistent with the 
fact that any single off-center vortex is unstable. Furthermore, a vortex state located at the center 
of the cloud first becomes locally stable as the rotational frequency increases. Finally our study 
implies the existence of hysteresis effects. 
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When rotated, a superfluid forms quantized vortex 
states. Numerous studies have examined in the past vor- 
tices in the traditional superfluid liquid Helium IV. Re- 
cently in some remarkable experiments vortices have also 
been created and observed in vapors of trapped ultracold 
atoms US 0,0. 

One of the basic questions in atomic systems is the for- 
mation and stability of vortex states. Because of the con- 
finement (which is typically harmonic) , there is a number 
of differences as compared to homogeneous superfluids. 
For example, the energy spectrum is discreet, the density 
is inhomogeneous and finally for harmonic confinement 
the frequency of rotation of the gas is limited by the trap 
frequency. 

Many studies have examined the physics of vortices in 
tr apped gases Hi,0 S 111313 El 111 13 13 E3 The 

basic picture is that above a critical frequency of rotation 
of the trap, one vortex state forms in the gas, while as 
the rotational frequency increases further, more vortices 
enter the cloud, eventually forming an array. One of the 
most fundamental and important questions is thus the 
way that the first vortex state forms. 

In the present study we develop a method which allows 
us to examine this problem in the limit of weak interac- 
tions, where the typical atom-atom interaction energy is 
smaller than the oscillator quantum of energy and one 
can restrict himself to the subspace of states in the low- 
est Landau level. In this limit the gas has a very peculiar 
property when the angular momentum per atom ranges 
between zero and unity, as the interaction e nerg y scales 
linearly with the angular momentum 0, UtI llStTTsl l2fj| , 
just like as in an ideal gas. This is an exact result within 
the lowest Landau level subspace of states. 

The linear behavior of the spectrum has important im- 
plications on the rotational properties of the gas, as at 
a critical frequency of rotation which is smaller than the 
trap frequency by a (small) amount that is proportional 
to the ratio between the interaction energy and the os- 
cillator energy, the gas is predicted to undergo a discon- 
tinuous phase transition from a non-rotating state to a 
state with one vortex at the center of the trap EH • 

However, as we show below, one needs to be careful 
with the effect of the truncation to the lowest Landau 
level wavef unctions. In the present study, using second- 
order perturbation theory we calculate the energy of the 



system to second order, where the energy no longer in- 
creases linearly with the angular momentum. From the 
derived dispersion relation we then examine the energy 
in the rotating frame, finding that the gas exhibits a non 
trivial and interesting behavior as function of the rota- 
tional frequency of the trap f2. 

Our study suggests that any single off-center vortex 
state is unstable. Furthermore, a vortex state that is 
located at the center of the cloud first becomes locally 
stable. We also predict that the gas comes to rest at a 
smaller than the one where rotation sets in (see Fig. 1), 
and the system exhibits hysteresis EH > which is a general 
characteristic of first-order phase transitions. The graphs 
in Fig. 1 show schematically the energy of the gas in the 
rotating frame for certain physically-relevant rotational 
frequencies which we calculate below. 

In our model we consider atoms interacting via a short- 
range effective interaction, 



V5 n t = Ub2*(ri-r i )/2. 



(1) 



Here Uq = A-kTi a/M is the strength of the effective two- 
body interaction, where a is the scattering length for elas- 
tic atom-atom collisions and M is the atom mass. We 
also consider a harmonic trapping potential of the form 

Muo 2 {p 2 + Az 2 )/2, 



V(p,z) 



(2) 



where p and z are cylindrical polar coordinates, u> is the 
oscillator frequency, and A is a dimensionless constant 
(h = M = ui = 1 from now on.) 

We consider weak interactions and strong confinement 
along the z axis, which is taken to be the axis of rotation. 
For weak interactions the corresponding dimensionless 
quantity is 7 = Na/d z which is our expansion parameter 
as it is assumed to be much smaller than unity. Here 
N is the total number of atoms and d z is the oscillator 
length along the axis of rotation. In this limit one can 
work within the subspace of the nodeless (n = 0) eigen- 
functions of the two dimensional harmonic oscillator, 



1pn,m(p, 4>) = 



7r(n + |jto|)! 



p |m| e™^L™(p 2 )e- p2 / 2 , (3) 



where n is the number of radial excitations, m is the 
quantum number of the angular momentum, and are 
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FIG. 1: Schematic diagram which shows the energy of the gas 
in the rotating frame as function of the angular momentum 
per atom Z for different frequencies of rotation Q. On the 
top graph, for Q — Q\, the derivative of the energy at Z = 1 
vanishes. For £1 = SI2 the energy at / = equals that at 
1 = 1. Finally, for £1 = Q3 the derivative of the energy at I — 
vanishes. Our study shows that fii < O2 < ^3, and also that 
the energy bents downward at both Z = and 1 = 1. All these 
observations have crucial consequences on the formation and 
stability of a single vortex state. 



the associated Laguerre polynomials. More specifically, 
the basis states are ip ,m(p, 4>) = p\ m \e im ' t 'e~ p2 / 2 / y/-n\m\\. 

The strong confinement along the axis of rotation, 
A ^> 1, implies that the cloud is in its lowest state of 
motion along this axis, and the problem thus becomes 
effectively two dimensional, as the degrees of freedom 
along the z axis are frozen out. Therefore, the order 
parameter ^(r) can be expanded in the product states 
*o,m(r) = ifo,m(p,<t>)<h( z )> where (po(z) is the ground 
state of the one dimensional harmonic oscillator. 

Under the above conditions, as shown initially by 
Bertsch and Papenbrock [17| using numerical diagonal- 
ization, the interaction energy of the gas in the lowest 
state is, for 2 < L < N, 
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FIG. 2: The five diagrams which contribute to the energy of 
the gas to order 7 2 and up to Z 3/2 , for 1 < L < N. The 
unprimed (primed) numbers denote the angular momentum 
m of the states with zero (any nonzero) radial excitations, 
n = (n > 1). For the states with m < 0, n can also be zero. 
The dotted lines denote the interaction. 



where L is the total angular momentum, and therefore it 
varies linearly with L. References have shown ana- 
lytically that this equation is exact to first order in 7. 

Therefore, if F = E — Ltt is the total energy of the 
gas in the rotating frame, where E = L + £l,n is the 
total energy in the rest frame then to first order in 7 
(henceforth the energy is measured with respect to the 
energy of the lowest state and terms of order l/N are 
neglected), F/N = (1 - 7/2 v^- Q) I, where I = L/N. 
The above equation implies that the critical frequency 
for rotation is 



f}W = 1 - 7 /2V2^. 



(5) 



At this value of £1 Butts and Rokhsar 01 predict that 
the gas undergoes a discontinuous transition from a non 
rotating state to a state with a vortex located at the 
center of the trap. In the present study we calculate 
the energy to next order for < I < 1 examining how 
this picture is modified. On the other hand, for I > 1 
where more than one vortices are present, calculation of 
the energy to lowest order in 7 suffices to determine their 
stability. 

We use now second-order perturbation theory to cal- 
culate the interaction energy to next order, i.e., to 7 2 . 
In a similar method Ref. [2(| has calculated the low-lying 
excitations of the system in the limit of low angular mo- 
mentum. In our study we use the results of Ref. [Is| 
which has studied both regimes of low (1 -c L <C N) and 
high (1 <C N — L -C N) angular momentum. Starting 
with 1 « L < JV, since we need both the slope, as well 
as the curvature of the energy (which is given by a term 
of order I 3 ' 2 in this case), the order parameter is, up to 
the desired order. 



* = c *o.o + c 2 *o,2 + c 3 *o„ 



(6) 



where |c | 2 = 1-1/2 + Z 3 / 2 /3, |c 2 | 2 = 1/2 - Z 3 / 2 , and 
I C3 1 2 = 2Z 3 / 2 /3. The five diagrams which contribute to 
the energy to second order in 7 and up to Z 3 / 2 are shown 
in Fig. 2. For example, the contribution of the top left 
diagram to the energy per particle is equal to 

N) |<*o,o,*o,o|M„t|*o,o,^o)| 2 ^ l4bT^o 2 



£l,n = l[N - L/2 - l)/y/to, 



(4) 
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FIG. 3: The five diagrams which contribute to the energy of 
the gas to order -f 2 and up to I 2 , for 1 <C N — L <C N. 
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FIG. 4: The nine diagrams which contribute to the energy of 
the gas to order y 2 and to I 2 , for 1 < N - L < N. 



{l-l/2 + l 3/2 /3) 3 (NU ) 2 J2 



\I 



0,ni2 
0,0 I 



£n,0 — £0,0 



-,(7) 



where I 'q = J ^o.o^o.o^o.o^n.o dr is the overlap in- 
tegral between the corresponding states. This always 
involves three states with zero radial excitations and 
one state with n radial excitation(s). Also t n>m = 
2n + \m\ + 1 + A/2 is the eigenenergy of the states 
*n,m( r ) = ^n,m(p>0)0o(^) and N a = A^|c | 2 is the oc- 
cupancy of the state ^o.o- In Eq. Q there is a factor of 
1/2 from the interaction [Eq. JTJ] that is cancelled by a 
factor of 2 that comes from the symmetrization of the 
wavefunction. 

The sum in Eq. J7J is over all the excited states with 
m = and n > 1 (in other cases where intermediate 
states with negative values of m are involved, n > 0.) 
However, the series converges rapidly because of the over- 
lap integrals which decrease with increasing n. Consid- 
ering the first fifteen excited states, n = 1, . . . , 15 we find 
that 

E/N m (l- j/2\/2nj I 

- (0.0916 - 0.0883 1 + 0.0156 Z 3/2 ) 7 2 . (8) 

As a final step we express E/N as function of the ex- 
pectation value of the angular momentum per parti- 
cle (I) = (*|£|*)/(*|§). Here * = E m c m * m , with 
* TO = *o,™ + 7E„/o d n ,m^n,m being the perturbed ba- 



sis states. The corrections of order 7 in the basis states 
introduce a correction of order 7 2 in (I) — I. 

To calculate (I) we start from the Gross- 
Pitaevskii equation and make use of the or- 
thogonality between the states ^ n ,m finding 
dn, m = -271-/ |* ,m| 2 *o,m*;, m rfr/7i. This for- 
mula implies that (I) = l[l - 0.0405 7 2 - 0.0018 7 2 1 1 / 2 ], 
or I = (Z)[l + 0.0405 7 2 + 0.0018 7 2 (l) 1/2 }. Combining 
this result with Eq. iJSJ we find 

E/N » (l - 7 /2V2^) (I) 

- (0.0916 - 0.1288 (I) + 0.0138 (Z) 3/2 ) 7 2 . (9) 

From Eq. |J3J we conclude that the critical frequency of 
rotation is, to order 7 2 , 

n ( V = 3 = 1 - 7 /2V2^ + 0.1288 7 2 + 0(7 3 ). (10) 

Also E/N (as well as F/N) has a downward curvature. 

Following a similar method we perform the same anal- 
ysis for 1 -c N — L -C N. There is no term of order I 3 / 2 
in this case and the two leading terms are of order I and 
I 2 , where 1 = 1 — 1. The corresponding order parameter 

is m 



0,0 



Cl*0,l+C 2 * ,2 + C 3 * 



0,3, 



(11) 



c 2 



where |c | 2 =21- 3Z 2 /2, | Cl | 2 = 1-31 + 271 2 
1-9P/4, and |c 3 | 2 =3f 2 /8. 

The five diagrams which contribute to the interaction 
energy linearly and quadratically in I are shown in Fig. 3, 
while the ones in Fig. 4 contribute only quadratically. 
Again, considering the first fifteen excited states we find 



E/N w ( 1 - 7/2V2vrJ (1 - I) 

' 0.0111 + 0. 1024 Z + 0. 7654 f 2 ) 7 2 . (12) 



In this case (Z) = l[l + 0.0780 7 2 - 0.0811 1-/ 2 }, or I = 
(Z)[l - 0.0780 7 2 + 0.0811 (I) 7 2 ]. The energy, expressed 
as function of (l) is therefore 



E/N re (1 - 7 /2V27t) (1-(Z) 

- 1 0.0111 + 0.0244 (1) + 0.8465 (Z~) 2 J 7 2 . (13) 

The frequency at which the derivative vanishes is thus 

Oi = 1 - 7 /2V2^ + 0.0244 7 2 + 0( 7 3 ). (14) 

Again, E/N (and F/N) has a downward curvature. 

Finally, another relevant frequency of rotation is the 
one where the energy in the rotating frame for (I) = 
equals that for (I) = 1, and this turns out to be 

n 2 = 1 -7/2V2t7+ 0.0805 7 2 + 0(7 3 ), (15) 
in agreement with Ref. [2l| . 
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Having calculated the energy as function of the angu- 
lar momentum in the two limiting cases, we can extract 
valuable information about the formation and stability of 
a single vortex state. Since our Hamiltonian is rotation- 
ally invariant, it commutes with the angular momentum 
and therefore the angular momentum is a good quantum 
number. In that respect, any configuration is stable, how- 
ever the interesting question is the stability against weak 
perturbations. 

For example, in the presence of a small thermal com- 
ponent in the gas which interacts with the condensate 
exchanging angular momentum and energy with it, a vor- 
tex state is stable/metastable as long as its energy in the 
rotating frame has an absolute/local minimum [22j |. In 
the limit of weak interactions that we consider here, the 
energy is dominated by the oscillator energy and the en- 
ergy barriers are small. As a result, these systems cannot 
support persistent currents [i.e., F(l) does not have any 
metastable minimum at any I ^ when Q = 0]. Still, 
when Q ^ 0, F(l) develops in principle local minima 
and in what follows, we consider a vortex state as a sta- 
ble configuration provided that F(l) has a local/absolute 
minimum (henceforth I is to be identified as (/)). 

The angular momentum plotted on the horizontal axis 
of Fig. 1 is directly related to the position of the vortex, 
since for an off-center vortex state < I < 1 [23j . There- 
fore, as I increases in Fig. 1 the vortex moves from an 
infinite distance away from the trap (I = 0) to its center 
(1 = 1). 

Our results here are exact for weak interactions and for 
values of the angular momentum close to zero and unity. 
Still, our expansion strongly suggests that the schematic 
form of F(l) shown in Fig. 1 extends over all the inter- 
mediate values of the angular momentum, without any 
local minima in between. 

As shown in Fig. 1 for any frequency f2 between fii and 



0,2, F(l) has a local minimum at I = 1. In other words 
a vortex that is located at the center of the trap first 
becomes stable locally with increasing f2. This behavior 
is qualitatively the same as in Ref. ^(| where the opposite 
limit of strong interaction was considered. 

According to our study, if one first cools down be- 
low the condensation temperature and then rotates, the 
gas will undergo a discontinuous transition from a non- 
rotating state to a state with a vortex in the middle of the 
cloud for r2 — r?3 . In the reverse process as f2 decreases, 
the gas will make a discontinuous transition from / = 1 
to / = at an = Qj, with Sli < fl^. 

In the other physically-relevant situation, if one first 
rotates with some f2 = Qq and then cools down below the 
condensation temperature, the gas may actually reach 
the state with a (stable) centered vortex for any value of 
S7o between fii and f^. 

It is crucial to mention that the slope of F(l) increases 
discontinuously as one crosses the point I = 1 by an 
amount of order 7, equal to ll7/32v / 27r [18j (not shown 
in Fig. 1). This fact guarantees that F(l) increases for 
I > 1 for all the values of fl considered above. 

To summarize, the spectrum of a weakly-interacting 
Bose-Einstein condensate that is confined in a rotating 
harmonic trap implies that a vortex state located at the 
center of the cloud becomes first locally stable. On the 
other hand, any single off-center vortex is unstable. Fi- 
nally we predict hysteresis as the frequency of rotation 
of the trap varies. 
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